• Computational Bayesian methods

• Spatio-temporal statistics

• Science-motivated statistical modeling

1. Spatio-temporal reconstruction of climate from pollen
   1. Non-linear functional relationship
      • Computationally challenging
      • Non-linear and non-Gaussian relationships among data
   2. Using data/parameter augmentation to improve computation
2. Recursive Bayesian inference for latent spatial data
   • Reducing computational cost through parallelization
   • Better Markov Chain Monte Carlo (MCMC) mixing

• Globally distributed.

• Trees are long-lived – little evolutionary drift during the Holocene.

• Current species distributions with respect to climate are well understood.

• Many databases of pollen data.
  • Neotoma.
  • European Pollen Database.
  • Global Pollen Database.
  • NOAA Pollen Database.

• Due to challenging data likelihood.

• Extremely noisy relationship between data and process.

  ---

  • No obvious patterns when plotting the data.

  ---

  • Inverse prediction problem.

• Lack of good statistical software.

  ---

  • gitHub packages pollenReconstruction, pgR, and BayesMRA in development.

• Link the observed pollen counts to climate states during the modern period.

• Use the learned relationship to predict unobserved climate state.

• Generate climate histories that are local to the site of interest with uncertainty.

• Spatially explicit reconstructions of climate variables is important.

  ---

  • Many important ecological questions are local.

  ---

  • Predictions at all locations remove the effects of sampling bias in paleoclimate reconstruction.

• Prior work – 4 sites and total compute time of approximately 28 hours.

  ---

  • Currently: 363 sites with
    • non-linear response: compute time \(\approx\) 48 hours.
    • Pólya-gamma linear response: compute time less than 2 hours.

• Sediment samples from a lake.

• Take 1cm\(^3\) cubes along the length of the sediment core.

• In each cube, researcher counts the first \(M\) pollen grains and identifies to species/OTU.

• Raw data are counts of each species.

• Compositional count data.

For the \(i\)th observation at location \(\mathbf{s}\) and time \(t\),

\[ \begin{align*} \mathbf{y}_i \left( \mathbf{s}, t \right) & = \left( y_{i1} \left( \mathbf{s}, t \right), \ldots, y_{id} \left( \mathbf{s}, t \right) \right)' \end{align*} \]

is a \(d\)-dimensional compositional count observation.

• \(y_{ij} \left( \mathbf{s}, t \right)\) is the count of species \(j\) in the \(i\)th sample at location \(\mathbf{s}\) and time \(t\).

\[ \begin{align*} \mathbf{y}_i \left( \mathbf{s}, t \right) | \mathbf{p}\left( \mathbf{s}, t \right) & \sim \operatorname{Multinomial} \left( M_i \left( \mathbf{s}, t \right), \mathbf{p}\left( \mathbf{s}, t \right) \right) \end{align*} \]

• \(M_i \left( \mathbf{s}, t \right) = \sum_{j=1}^d y_{ij}\left( \mathbf{s}, t \right)\) is the total count observed (fixed and known) for observation \(i\) at location \(\mathbf{s}\) and time \(t\).

• Compositional count data.

  ---

  • Total count is not informative of the absolute composition.

  ---

  • Observation informative of the relative proportions \(p_{j} \left( \mathbf{s}, t \right)\) only.

• The pollen data are highly variable and overdispersed.

  ---

• Mixture over a Dirichlet distribution.

  ---

\[ \begin{align*} \mathbf{p}\left( \mathbf{s}, t \right) | \boldsymbol{\alpha}\left( \mathbf{s}, t \right) & \sim \operatorname{Dirichlet} \left( \boldsymbol{\alpha}\left( \mathbf{s}, t \right) \right) \end{align*} \]

• Marginalize out \(\mathbf{p} \left( \mathbf{s}, t \right)\).

\[ \begin{align*} \mathbf{y}_i\left( \mathbf{s}, t \right) | \boldsymbol{\alpha}\left( \mathbf{s}, t \right) & \sim \operatorname{Dirichlet-Multinomial} \left( M_i\left( \mathbf{s}, t \right), \boldsymbol{\alpha}\left( \mathbf{s}, t \right) \right) \end{align*} \]

• Model the Dirichlet-multinomial random effect using the log link function:

\[ \begin{align*} \operatorname{log} \left( \boldsymbol{\alpha} \left( \mathbf{s}, t \right) \right) & = \mathbf{z}\left( \mathbf{s}, t \right) \boldsymbol{\beta}. \end{align*} \]

• \(\mathbf{z}\left( \mathbf{s}, t \right)\)’ is a \(q\)-dimensional vector of climate variables.

• \(\boldsymbol{\beta}\) is a \(q \times d\) dimensional matrix of regression coefficients.

• Vegetation response to climate is non-linear.

\[ \begin{align*} \operatorname{log} \left( \boldsymbol{\alpha} \left( \mathbf{s}, t \right) \right) & = \mathbf{B} \left( \mathbf{z}\left( \mathbf{s}, t \right) \right) \boldsymbol{\beta} \end{align*} \]

• \(\mathbf{B} \left( \mathbf{z}\left( \mathbf{s}, t \right) \right)\) is a basis expansion of the covariates \(\mathbf{z}\left( \mathbf{s}, t \right)\).
  • Use B-splines or Gaussian Processes as a basis.
  • \(\mathbf{z}\left( \mathbf{s}, t \right)\) is unknown for \(t \neq 1\).
    
  • Computationally challenging.

Nolan C, Tipton JR, Booth RK, Hooten MB, Jackson ST (In Press)., “Comparing and improving methods for reconstructing peatland, water-table depth from testate amoebae.” The Holocene.

Hefley TJ, Broms KM, Brost BM, Buderman FE, Kay SL, Scharf HR, Tipton, JR, Williams PJ, Hooten MB (2017). “The basis function approach for, modeling autocorrelation in ecological data.” Ecology, 98(3),, 632-646.

• \(\mathbf{z} \left( \mathbf{s}, t \right)\)s are observed only at \(t\) = 1.

• Calibration: Estimate \(\boldsymbol{\beta}\) using: \[ \begin{align*} \mathbf{y} \left(1\right) & = \left( \mathbf{y} \left( \mathbf{s}_1, 1 \right), \ldots, \mathbf{y} \left( \mathbf{s}_n, 1 \right) \right)' \\ \mathbf{z} \left(1\right) & = \left( \mathbf{z} \left( \mathbf{s}_1, 1 \right), \ldots, \mathbf{z} \left( \mathbf{s}_n, 1 \right) \right)' \end{align*} \]

• Reconstruction:
  • Use estimated \(\boldsymbol{\beta}\)s and fossil pollen \(\mathbf{y} \left( t \right)\) to predict unobserved \(\mathbf{z}\left( t \right)\).

• Estimate the latent process \(\mathbf{z} \left( \mathbf{s}, t \right)\).

• Options:
  1. continuous vs. discrete space (geostatistical vs. CAR models).
  2. continuous vs. discrete time (stochastic process vs. AR models).

• For now, we focus on continuous space and discrete time.

• For \(\mathbf{z} \left(t \right) = \left( \mathbf{z} \left(\mathbf{s}_1, t \right)', \ldots, \mathbf{z} \left(\mathbf{s}_n, t \right)' \right)\), we assume:

\[ \begin{align*} \mathbf{z} \left(t \right) - \mathbf{X} \left( t \right) \boldsymbol{\gamma} & = \mathbf{M}\left(t\right) \left( \mathbf{z} \left(t-1 \right) - \mathbf{X} \left( t \right) \boldsymbol{\gamma} \right) + \boldsymbol{\eta} \left(t \right) \end{align*} \]

• \(\mathbf{M}(t) = \rho \mathbf{I}_n\) is a propagator matrix.

• \(\mathbf{X} \left(t \right) \boldsymbol{\gamma}\) are the fixed effects from covariates like latitude, elevation, etc.

• \(\boldsymbol{\eta} \left( t \right) \sim \operatorname{N} \left( \mathbf{0}, \mathbf{R}\left( \boldsymbol{\theta} \right) \right)\).

• \(\mathbf{R} \left( \boldsymbol{\theta} \right)\) is a Mátern spatial covariance matrix with parameters \(\boldsymbol{\theta}\).

• Define a set of spatial knot locations \(\mathbf{s}^{\star} = \left\{ \mathbf{s}_1^{\star}, \ldots, \mathbf{s}_m^{\star} \right\}\).

\[ \boldsymbol{\eta}^{\star} \left( t \right) \sim \operatorname{N} \left( \mathbf{0}, \mathbf{R}^{\star}\left( \boldsymbol{\theta} \right) \right) \]

• \(\mathbf{R}^{\star}\left( \boldsymbol{\theta} \right)\) is the spatial covariance defined at the knot locations \(\mathbf{s}^{\star}\).

• Linear interpolator from observed locations \(\mathbf{s}\) to knot locations \(\mathbf{s}_j^{\star}\) is \[ \mathbf{r} \left(\mathbf{s}, \mathbf{s}_j^{\star} \right) \mathbf{R}^{\star}\left( \boldsymbol{\theta} \right)^{-1} \] where \(\mathbf{r} \left(\mathbf{s}, \mathbf{s}_j^{\star} \right) = \operatorname{Cov} \left(\mathbf{s}, \mathbf{s}_j^{\star} \right)\)

Banerjee S, Gelfand AE, Finley AO, Sang H (2008). “Gaussian predictive, process models for large spatial data sets.” Journal of the Royal, Statistical Society: Series B (Statistical Methodology), 70(4),, 825-848.

• \(\boldsymbol{\eta} \left( t \right) \approx \mathbf{r} \left(\mathbf{s}, \mathbf{s}^{\star} \right) \mathbf{R}^{\star}\left( \boldsymbol{\theta} \right)^{-1} \tilde{\boldsymbol{\eta}} \left( t \right)\).

• The predictive process can be shown to be the optimal predictor of the parent process \(\boldsymbol{\eta} \left( t \right)\) of dimension \(m\)

• The dynamic climate process becomes

\[ \begin{align*} \mathbf{z} \left(t \right) - \mathbf{X} \left( t \right) \boldsymbol{\gamma} & \approx \mathbf{M}\left(t\right) \left( \mathbf{z} \left(t-1 \right) - \mathbf{X} \left( t \right) \boldsymbol{\gamma} \right) + \mathbf{r} \left(\mathbf{s}, \mathbf{s}^{\star} \right) \mathbf{R}^{\star}\left( \boldsymbol{\theta} \right)^{-1} \boldsymbol{\tilde{\eta}} \left(t \right) \end{align*} \]

• Each fossil pollen observation includes estimates of time uncertainty.
  • The time of the observation is uncertain.
  • Weight the likelihoods according to age-depth model.
  • Posterior distribution of ages.

• For each observation fossil pollen observation an age-depth model gives a posterior distribution over dates.
  • Define \(\omega \left(\mathbf{s}, t \right)\) as P(age \(\in (t-1, t)\)).
  • \([\mathbf{y} \left( \mathbf{s}, t \right) | \boldsymbol{\alpha} \left( \mathbf{s}, t \right) ] = \prod_{t=1}^T [\mathbf{y} \left( \mathbf{s}, t \right) | \boldsymbol{\alpha} \left( \mathbf{s}, t \right)]^{\omega_\left(\mathbf{s}, t \right)}\).

• Non-spatial code available as a gitHub package.

devtools::install_github("jtipton25/BayesComposition")

• Includes options for multiple models including:
  • mixture models.
  • different likelihoods and link functions.
  • correlations in functional response.

• Code in C++ using Rcpp package for computation speed.

1. Fit the calibration model.

2. Use posterior distribution from stage (1) to generate predictions of climate independent in space and time.

3. Smooth the posterior distribution from stage (2) using dynamic linear model.

   • Use posterior mean estimates which does not fully quantify model uncertainty.

   • Goal: Use recursive Bayesian ideas from end of talk.

Hooten MB, Johnson DS, Brost BM (2019). “Making Recursive Bayesian, Inference Accessible.” The American Statistician, 1-10.

• Single stage Bayesian model.

• High-dimensional spatio-temporal process.
  • Inefficient to sample with block Metropolis.
  • Poor mixing of MCMC chains.

• Non-linear transformation in the data model.
  • Difficult to use Kalman Filtering.

• Particle Filtering Methods.
  • Difficult to implement, suffer from degeneracy.
  • Posterior can be multi-modal.

• Elliptical Slice Sampling.
  • Assumes a Gaussian prior.
  • Difficult to implement in a DLM framework.
  • Requires no tuning.
  • Efficiently samples in high dimensions.
  • Easily explores multiple modes.

• Adaptive block Metropolis within Gibbs and Elliptical Slice Sampling algorithms.

• Highly multi-modal posterior is efficiently explored within the sampler.

Murray I, Adams RP, MacKay DJC (2010). “Elliptical slice sampling.” In, AISTATS, volume 13, 541-548.

• For Gaussian process expansion of \(\mathbf{z} \left( \mathbf{s}, t \right)\):
  • The latent climate states are inputs into the covariance function.
  • Covariance input locations (and distance) are random.
  • Unique computational challenge.
    • Total computational cost is prohibitive \(O(d \frac{(Tn)^3}{3})\).
• Proposed solution:
  • Predictive process representation.
  • Reduced computation cost \(O(d T \frac{m^3}{3})\)

Tipton JR, Hooten MB, Nolan C, Booth RK, McLachlan J (2019)., “Predicting paleoclimate from compositional data using multivariate, Gaussian process inverse prediction.” The Annals of Applied, Statistics, 13(4), 2363-2388.

\[\begin{align*} \mathbf{y}(\mathbf{s}, t) & \sim \operatorname{Multinomial}(M(\mathbf{s}, t), \pi_{SB}(\boldsymbol{\eta}(\mathbf{s}, t))) \end{align*}\]

• \(\mathbf{y}(\mathbf{s}, t)\) is a \(J\) dimensional vector of counts at site \(\mathbf{s} \in \mathcal{D}\) and time \(t \in \{1, \ldots, n_t\}\)

• \(M(\mathbf{s}, t) = \sum_{j=1}^J \mathbf{y}_j(\mathbf{s}, t)\) is the total count observed.

• Link the underlying climate states to the probability of observing species \(j\) through the latent variable \(\eta_j(\mathbf{s}, t)\)

• \(\pi_{SB}(\boldsymbol{\eta}(\mathbf{s}, t))\) is a stick breaking transformation.

  ---

  • The map \(\pi_{SB}: \mathcal{R}^{J-1} \rightarrow \Delta^{J}\) transforms the \(J-1\) dimensional vector \(\boldsymbol{\eta}(\mathbf{s}, t)\) to the \(J\) dimensional unit simpled \(\mathcal{\Delta}^{J}\).

• Other maps to the unit simplex could be used (i.e., multi-logit), but the stick-breaking map will be important later.

\[\begin{align*} [\mathbf{y} | M, \pi_{SB}(\boldsymbol{\eta})] & =\operatorname{Multinomial}(\mathbf{y} | M, \pi_{SB}(\boldsymbol{\eta})) \\ & = \prod_{j=1}^{J-1} \operatorname{binomial}(y_j | M_j, \tilde{\pi}_j) \\ & = \prod_{j=}^{J-1} {M_j \choose y_j} \frac{(e^{\eta_j})^{y_j}}{(1 + e^{\eta_j})^{M_j} } \end{align*}\]

• \(M_j = M - \sum_{k < j} M_k\).

• Define the partial probabilities \(\tilde{\pi}_j = \pi_{SB}(\boldsymbol{\eta})_j\) using the stick-breaking representation.

\[\begin{align*} \frac{(e^{\eta_j})^{y_j}}{(1 + e^{\eta_j})^{M_j} } & = 2^{-M_j} e^{\kappa(y_j) \eta_j} \int_0^\infty e^{- \omega_j \eta_j^2 / 2 } [\omega_j | M_j, 0] \,d \omega \end{align*}\]

• The integral identity is proportional to the product representation of the multinomial distribution

• The density \([\omega_j | M_j, 0]\) is a Pólya-gamma distribution \(PG(b, c)\)

• With a prior \([\eta_j]\), joint density \([y_j, \eta_j]\) is

\[\begin{align*} [\eta_j, y_j] & = [\eta_j] {M_j \choose y_j} \frac{(e^{\eta_j})^{y_j}}{(1 + e^{\eta_j})^{M_j} }\\ & = \int_0^\infty [\eta_j] {M_j \choose y_j} 2^{-M_j} e^{\kappa(y_j) \eta_j} e^{- \omega_j \eta_j^2 / 2 } [\omega_j | M_j, 0] \,d \omega \end{align*}\]

where the integral defines a joint density over \([\eta_j, y_j, \omega_j]\).

• Using this integral representation, we have

  \[\begin{align*} \omega_j | \eta_j, y_j & \sim \operatorname{PG( M_j, \eta_j)} \end{align*}\]

  which can be sampled using the exponential tilting property of the Pólya-gamma distribution

• If \([\eta_j]\) is Gaussian, then \([\eta_j | \omega_j, y_j]\) is also Gaussian which enables conjugate sampling.

• There is a cost:

  ---

  • Requires sampling the \(\omega_j\)s and these are computationally expensive.

  ---

  • However, these are also embarrassingly parallel.

  ---

  • Efficiently sampled using openMP parallelization.

• For all examples tested so far with reduced-rank spatial processes, sampling Pólya-gamma random variables is the limiting computational cost.

\[\begin{align*} \eta_j(\mathbf{s}, t) = \color{blue}{\beta_{0j}} + \color{blue}{\beta_{1j}} \left( \mathbf{x}'(\mathbf{s}, t) \boldsymbol{\gamma} + \mathbf{w}'(\mathbf{s}) \boldsymbol{\alpha}(t) \right) + \color{blue}{\varepsilon(\mathbf{s}, t)} \end{align*}\]

• \(\color{blue}{\beta_{0j}}\) and \(\color{blue}{\beta_{1j}}\) are regression coefficients with respect to the climate state \(\mathbf{Z}(\mathbf{s}, t) = \left( \mathbf{x}'(\mathbf{s}, t) \boldsymbol{\gamma} + \mathbf{w}'(\mathbf{s}) \boldsymbol{\alpha}(t) \right)\).

• \(\color{blue}{\varepsilon(\mathbf{s}, t)} \stackrel{iid}{\sim} \operatorname{N}(0, \sigma^2_j)\) models overdispersion relative to the linear regression response.

\[\begin{align*} \eta_j(\mathbf{s}, t) = \beta_{0j} + \beta_{1j} \left( \color{red}{\mathbf{x}'(\mathbf{s}, t) \boldsymbol{\gamma}} + \color{purple}{\mathbf{w}'(\mathbf{s}) \boldsymbol{\alpha}(t)} \right) + \varepsilon(\mathbf{s}, t) \end{align*}\]

\(\color{red}{\mathbf{X}(t) = \begin{pmatrix} \mathbf{x}'(\mathbf{s}_1, t) \\ \vdots \\ \mathbf{x}'(\mathbf{s}_n, t) \end{pmatrix}}\) are fixed covariates (elevation, latitude, etc.).

• We assume \(\color{red}{\mathbf{X}(t) \equiv \mathbf{X}}\) for all \(t\) although temporally varying covariates are possible (volcanic forcings, Milankovitch cycles, etc.).

• \(\color{purple}{\mathbf{W} = \begin{pmatrix} \mathbf{w}'(\mathbf{s}_1) \\ \vdots \\ \mathbf{w}'(\mathbf{s}_n) \end{pmatrix}}\) are spatial basis functions with temporal random effects \(\color{purple}{\boldsymbol{\alpha}(t)}\).

• \(\mathbf{Z}_0 \sim \operatorname{N} (\color{red}{\mathbf{X}'(1) \boldsymbol{\gamma}} + \color{purple}{\mathbf{W} \boldsymbol{\alpha}(1)}, \sigma^2_0 \mathbf{I})\) is the observed modern climate state.

• Challenge: scaling to 1000s of sites and 15,000 years (\(\approx 60\) time increments)

• Sparse, multiresolution representation with a dynamic linear model

\[\begin{align*} \color{purple}{\mathbf{w}'(\mathbf{s}) \boldsymbol{\alpha}(t)} = \color{purple}{\sum_{m=1}^M \mathbf{w}_m'(\mathbf{s}) \boldsymbol{\alpha}_m(t)} \end{align*}\]

• \(\color{purple}{\mathbf{w}_m'(\mathbf{s})}\) is a Wendland basis over resolution \(m\)

\[\begin{align*} \color{purple}{\boldsymbol{\alpha}_m(t)} & \sim \operatorname{N} \left( \mathbf{A}_m \color{purple}{\boldsymbol{\alpha}_m(t-1)}, \tau^2 \mathbf{Q}_m^{-1} \right) \end{align*}\]

• Assume the constraint \(\color{purple}{(\mathbf{W}_m \boldsymbol{\alpha}_m(t))'} \mathbf{1} = 0\) for all \(m\) and \(t\)

• \(\mathbf{A}_m\) is the evolution matrix (assume \(\mathbf{A}_m \equiv \rho \mathbf{I}\) for all \(m\))

• \(\mathbf{Q}_m\) is the precision matrix constructed from either a CAR or SAR process over the \(m\)th resolution

• Consider the Gaussian linear mixed model - first order (mean) representation:

\[ \begin{align*} \mathbf{y} & = \mathbf{X} \boldsymbol{\beta} + \mathbf{Z} \boldsymbol{\alpha} + \boldsymbol{\varepsilon} \end{align*} \]

• \(\mathbf{X}\) is a set of fixed effect covariates.

• \(\boldsymbol{\beta}\) are regression coefficients.

• \(\mathbf{Z}\) is a set of random effect covariates and/or basis functions.

• \(\boldsymbol{\alpha} \sim \operatorname{N}(\mathbf{0}, \boldsymbol{\Sigma}_{\alpha})\) are random effects.

• \(\boldsymbol{\varepsilon} \sim \operatorname{N}(\mathbf{0}, \sigma^2 \mathbf{I})\)

• First-order (mean) representation

\[ \begin{align*} \mathbf{y} & \sim \operatorname{N} \left(\mathbf{X} \boldsymbol{\beta} + \mathbf{Z} \boldsymbol{\alpha}, \sigma^2 \mathbf{I} \right) \\ \boldsymbol{\alpha} & \sim \operatorname{N} \left( \mathbf{0}, \boldsymbol{\Sigma}_{\alpha} \right) \\ \boldsymbol{\beta} & \sim \operatorname{N} \left( \mathbf{0}, \boldsymbol{\Sigma}_{\beta} \right) \end{align*} \]

• Second order (covariance) representation

\[ \begin{align*} \mathbf{y} & \sim \int \int \operatorname{N} \left(\mathbf{X} \boldsymbol{\beta} + \mathbf{Z} \boldsymbol{\alpha}, \sigma^2 \mathbf{I} \right) \operatorname{N} \left( \mathbf{0}, \boldsymbol{\Sigma}_{\alpha} \right) \operatorname{N} \left( \mathbf{0}, \boldsymbol{\Sigma}_{\beta} \right) \, \mathrm{d} \boldsymbol{\alpha} \, \mathrm{d} \boldsymbol{\beta} \\ & \sim \operatorname{N} \left(\mathbf{0}, \sigma^2 \mathbf{I} + \mathbf{Z} \boldsymbol{\Sigma}_{\alpha}\mathbf{Z'} + \mathbf{X} \boldsymbol{\Sigma}_\beta \mathbf{X}' \right) \end{align*} \]

• For Gaussian models, this integration is standard practice – improved mixing, convergence, and efficiency.

• For non-Gaussian models, brute force or approximations are used.

  • Sample the latent \(\boldsymbol{\alpha}\)s and \(\boldsymbol{\beta}\)s directly.

  • Integrated nested Laplace approximations (INLA).

  • Replace the non-Gaussian likelihood by an approximating distribution that is Gaussian – Variational Inference.

  • Congujate recursive Bayesian inference.

\[ \begin{align*} y \left( \mathbf{s} \right) & = \mathbf{X} \left( \mathbf{s} \right) \boldsymbol{\beta} + \eta \left( \mathbf{s} \right) + \varepsilon \left( \mathbf{s} \right) \\ \boldsymbol{\eta} & \sim \operatorname{N} \left( \mathbf{0}, \mathbf{C} \right) \end{align*} \]

• \(\boldsymbol{\eta} = \left( \eta \left( \mathbf{s}_1 \right), \ldots, \eta \left( \mathbf{s}_n \right) \right)'\).

• Spatial location index \(\mathbf{s}\).

• \(\mathbf{C}\left( \cdot, \cdot \right)\) is a spatial covariance function.

• \(\mathbf{C}\) is a spatial covariance matrix with \(ij\)th entry \(\mathbf{C}_{ij} = \mathbf{C}\left( \mathbf{s}_i, \mathbf{s}_j \right)\).

• Observations nearby are more correlated than observations further away – Toblers first law of geography.

• Let \(\mathbf{Z}\) be a set of spatial basis functions where \(\mathbf{Z} \mathbf{Z}' = \mathbf{C}\).

\[ \begin{align*} y \left( \mathbf{s} \right) & = \mathbf{X} \left( \mathbf{s} \right) \boldsymbol{\beta} + \eta \left( \mathbf{s} \right) + \varepsilon \left( \mathbf{s} \right) \\ & = \mathbf{X} \left( \mathbf{s} \right) \boldsymbol{\beta} + \mathbf{Z} \left( \mathbf{s} \right) \boldsymbol{\alpha} + \varepsilon \left( \mathbf{s} \right). \end{align*} \]

• \(\mathbf{Z}\) can be chosen to provide fast computation, reduced-rank computation, or other properties of interest if you use the approximation \(\mathbf{Z} \boldsymbol{\Sigma}_{\alpha} \mathbf{Z}' \approx \mathbf{C}\).

• Let \(\mathbf{C} = \mathbf{U} \mathbf{D} \mathbf{U}'\) be an eigen-decomposition. Then setting \(\mathbf{Z}= \mathbf{U}\) and \(\boldsymbol{\alpha} \sim \operatorname{N} \left( \mathbf{0}, \mathbf{D} \right)\) gives an equivalent model

\[ \begin{align*} \mathbf{y} & = \int \operatorname{N} \left( \mathbf{X} \boldsymbol{\beta} + \mathbf{U} \boldsymbol{\alpha}, \sigma^2 \mathbf{I} \right) \operatorname{N} \left( \mathbf{0}, \boldsymbol{D} \right) \, \mathrm{d} \boldsymbol{\alpha} \\ & = \operatorname{N} \left( \mathbf{X} \boldsymbol{\beta}, \sigma^2 \mathbf{I} + \mathbf{U} \mathbf{D} \mathbf{U}' \right) \\ & = \operatorname{N} \left( \mathbf{X} \boldsymbol{\beta}, \sigma^2 \mathbf{I} + \mathbf{C} \right) \\ \end{align*} \]

• Integrating out the latent mean structure improves model stability and results in better model fitting.

• However, integrating out the random effect (usually) requires a Gaussian likelihood.

• How can we perform integrated Bayesian inference to get the computational advantages of the second-order representation in non-Gaussian models?

Classical model

• Single stage (joint) model posterior

\[ \begin{align*} [\mathbf{z}, \boldsymbol{\theta}_z, \boldsymbol{\theta}_y | \mathbf{y}] & = \prod_{i=1}^N [y_i | z_i, \boldsymbol{\theta}_y] [\mathbf{z} | \boldsymbol{\theta}_z] [\boldsymbol{\theta}_y] [\boldsymbol{\theta}_z]. \end{align*} \]

• Specify full conditional distributions and sample.

• Metropolis-Hastings, Gibbs, slice sampling, Hamiltonian Monte Carlo, etc.

• \(\mathbf{z}\) is a spatial random process .

  • \(\mathbf{z} \sim \operatorname{N} (\mathbf{0}, \boldsymbol{C})\)

• Must sample the process model (random effect \(\mathbf{z}\)) directly – cannot integrate out the process variable.

• Computational complexity grows with \(O(n^3)\).

• Inefficient to sample from a complicated likelihood surface.

• Fit the model in two stages.

• Start with a simple first stage model.

• Use output of first stage model as input into second stage model.

• First stage model generates an importance sampling proposal distribution for the second stage.

• If properly chosen, can integrate out the latent process \(\mathbf{z}\) to reduce computational cost.

• First, fit a simple model:

\[ \begin{align*} [\tilde{z}_i | y_i] & = [y_i | \tilde{z}_i] [\tilde{z}_i ]. \end{align*} \]

• Assign each of the \(\tilde{z}_i\) independent transient priors \([\tilde{z}_i]\).

• Posterior should be overdispersed relative to the target posterior.

• Can choose \([\tilde{z}_i]\) to allow for analytic inference.

• Generate \(k_1 = 1, \ldots, K_1\) first stage samples \(\tilde{z}_i^{(k_1)}\).

1. Current state of variable \(z\) at the \(k\)th MCMC iteration is \(z^{(k)}\).

2. Propose \(z^{\star}\) from proposal distribution \([z^{\star}|z^{(k)}]\).

3. Calculate MH ratio \(a^{\star} = \frac{[y|z^{\star}] [z^{\star}] [z^{\star} | z^{(k-1)}]}{[y|z^{(k-1)}] [z^{(k-1)}] [z^{(k-1)} | z^{\star}]}\).

4. Set

\[ \begin{align*} z^{(k)} = \begin{cases} z^{\star} & \mbox{ with probability } \min(1, a^{\star}) \\ z^{(k-1)} & \mbox{otherwise}. \end{cases} \end{align*} \]

• Use the posterior samples \(\tilde{z}^{(k_1)}\) for \(k_1 = 1, \ldots, K_1\) from the first stage posterior distribution \([\tilde{z}_i | y_i]\) as a proposal for second stage estimation.

• MH update for the \(k_2\)th second stage iteration:

• Sample \(\tilde{z}_i^{\star}\) from the \(K_1\) first stage posterior samples from \([\tilde{z}_i | y_i ]\) (i.e., set \(\star\) to one of \(1, \ldots, K_1\)).

• MH-ratio is \[ \begin{align} \label{eq:MH} a_i^{(k_2)} & = \frac{\color{blue}{[y_i | \tilde{z}_i^{\star}]} [\tilde{z}_i^{\star} | \tilde{\mathbf{z}}_{-i}, \boldsymbol{\theta}_z^{(k_2-1)}] \color{red}{[\tilde{z}^{(k_2-1)}_i | y_i]} }{\color{red}{[y_i | \tilde{z}^{(k_2-1)}_i]} [\tilde{z}^{(k_2 - 1)}_i | \tilde{\mathbf{z}}_{-i}, \boldsymbol{\theta}^{(k_2-1)}] \color{blue}{[\tilde{z}^{\star}_i | y_i]} }. \end{align} \]

• Both \(\color{red}{[\tilde{z}^{(k_2-1)}_i | y_i]}\) and \(\color{blue}{[\tilde{z}^{\star}_i | y_i]}\) are samples from the first stage posterior:

• Replace these densities by their un-normalized first stage densities \(\color{red}{[y_i | \tilde{z}^{(k_2-1)}_i]} \color{orange}{[\tilde{z}^{(k_2-1)}_i]}\) and \(\color{blue}{[y_i | \tilde{z}^{\star}_i]} \color{orange}{[\tilde{z}_i^{\star}]}\).

• Then, the MH acceptance ratio becomes

\[ \begin{align} a_i^{(k_2)} & = \frac{\color{blue}{[y_i | \tilde{z}_i^{\star}]} [\tilde{z}_i^{\star} | \mathbf{z}_{-i}, \boldsymbol{\theta}_z^{(k_2-1)}] \color{red}{[y_i | \tilde{z}^{(k_2-1)}_i]} \color{orange}{[\tilde{z}^{(k_2-1)}_i]} } {\color{red}{[y_i | \tilde{z}^{(k_2-1)}_i]} [\tilde{z}^{(k_2 - 1)}_i | \tilde{\mathbf{z}}_{-i}, \boldsymbol{\theta}^{(k_2-1)}] \color{blue}{[y_i | \tilde{z}^{\star}_i]} \color{orange}{[\tilde{z}_i^{\star}]} }. \\ \nonumber & = \frac{[\tilde{z}_i^{\star} | \mathbf{z}_{-i}, \boldsymbol{\theta}_z^{(k_2-1)}] \color{orange}{[\tilde{z}^{(k_2-1)}_i]} }{[z^{(k_2 - 1)}_i | \mathbf{z}_{-i}, \boldsymbol{\theta}_z^{(k_2-1)}] \color{orange}{[\tilde{z}_i^{\star}]} }. \end{align} \]

• \(\color{orange}{[\tilde{z}_i^{\star}]}\) and \(\color{orange}{[\tilde{z}_i^{(k_2 - 1)}]}\) are the first-stage, transient prior distributions.

• The ratio \(\frac{\color{blue}{[y_i | z_i^{\star}]} \color{red}{[y_i | \tilde{z}^{(k_2-1)}_i]} }{ \color{red}{[y_i | z^{(k_2-1)}_i]} \color{blue}{[y_i | \tilde{z}^{\star}_i]} } = 1\) so this drops out of the MH ratio.

• Notice the following doesn’t use the data:

\[ \begin{align} a_i^{(k_2)} & = \frac{[z_i^{\star} | \mathbf{z}_{-i}, \boldsymbol{\theta}_z^{(k_2-1)}] [\tilde{z}^{(k_2-1)}_i] }{[z^{(k_2 - 1)}_i | \mathbf{z}_{-i}, \boldsymbol{\theta}^{(k_2-1)}] [\tilde{z}_i^{\star}] }. \end{align} \]

• If the data model \([y_i|z_i]\) is non-Gaussian but \([\mathbf{z}_i|\boldsymbol{\theta}_z]\) is Gaussian, we can integrate out the random effect.

• Allows for fast, efficient sampling.

• Negative – MH: many of the proposed values from \([\tilde{z}_i | y_i ]\) will not have a high probability of acceptance.

• Importance sample over all of the first stage samples \([\tilde{z}_i | y_i ]\).

• Non-normalized importance weights:

\[ \begin{align*} w_i^{(k_1, k_2)} & = [z_i^{(k_1)} | \mathbf{z}_{-i}, \boldsymbol{\theta}_z^{(k_2-1)}] [\tilde{z}^{(k_2-1)}_i]. \end{align*} \]

• Each of these weights requires evaluation of the density \([z_i^{(k_1)} | \mathbf{z}_{-i}, \boldsymbol{\theta}_z^{(k_2-1)}]\) which can be computationally costly and is evaluated \(N\) times per MCMC sample.

• Most of the first-stage samples are not likely in the second stage.

• Wasted computational effort to evaluate weights for unlikely samples.

• Subsampled importance samples:

  • Assign a latent variable \(\delta_i^{(k_1)}\) and only evaluate the importance weights when \(\delta_i^{(k_1)} = 1\).

  • Embarrassingly parallel update.

• Instead of using the first stage posterior samples from \([\tilde{z}(\mathbf{s}_i) | y(\mathbf{s})_i]\) as a MH proposal, apply a density estimation method (e.g., kernel density estimation, Dirichlet process estimation, etc.) and use the estimated density to generate proposals.

• This can be done relatively easily using either an accept-reject algorithm or an importance sampler.

• In the Bayesian spatial glm, the process model is of the form

\[ \begin{align*} y \left( \mathbf{s}_i \right) & = g^{-1} \left( z \left( \mathbf{s}_i \right) \right) \\ z \left( \mathbf{s}_i \right) & = \mathbf{X} \left( \mathbf{s}_i \right) \boldsymbol{\beta} + \eta \left( \mathbf{s}_i \right). \end{align*} \]

• \(g\left( \cdot \right)\) is an appropriately chosen link function.

• \(\mathbf{X} \left( \mathbf{s}_i \right)\) is a set of fixed effect covariates.

• \(\boldsymbol{\beta}\) are regression coefficients.

• \(\boldsymbol{\eta} = \left( \eta \left(\mathbf{s}_1 \right), \ldots, \eta \left( \mathbf{s}_N \right) \right)' \sim \operatorname{N} (\mathbf{0}, \mathbf{C})\) is a random spatial process.

• In general, sampling a high-dimensional latent effect like \(\boldsymbol{\eta}\) in BHMs is a challenging and unsolved problem.

• Typical solution for the Gaussian case – integrate out the latent variable and model the covariance.

• Not possible in non-Gaussian models (GLMs).

• In the two-stage model, the data distribution \([y(\mathbf{s})_i | z(\mathbf{s}_i), \boldsymbol{\theta}_y]\) does not appear in the Metropolis ratio .

  • In the second stage posterior, one can directly integrate out the latent random variable \(\boldsymbol{\eta}\) and there is no requirement to sample a high-dimensional random variable which improves MCMC mixing.

• Let \(\mathbf{s}_i\) represent a spatial location and \(y(\mathbf{s}_i)\) be the observation.

\[ \begin{align*} y(\mathbf{s}_i) | \lambda(\mathbf{s}_i) & \stackrel{iid}{\sim} \operatorname{Poisson}\left( \lambda \left( \mathbf{s}_i \right) \right), \end{align*} \]

where the latent intensity \(\boldsymbol{\lambda} = \left( \lambda \left( \mathbf{s}_1 \right), \ldots, \lambda \left(\mathbf{s}_n \right) \right)'\) is modeled using a log-link function with the distribution

\[ \begin{align*} \log(\boldsymbol{\lambda}) & = \mathbf{X} \boldsymbol{\beta} + \boldsymbol{\eta} \\ \boldsymbol{\eta} | \tau^2, \boldsymbol{\theta} & \sim \operatorname{N} \left( \mathbf{0}, \tau^2 \mathbf{R} \left( \boldsymbol{\theta} \right) \right), \end{align*} \] where \(\tau^2\) is a variance parameter and \(\boldsymbol{\theta}\) is a set of parameters for the spatial correlation matrix \(\mathbf{R} \left( \boldsymbol{\theta} \right)\).

• Note: \(\boldsymbol{\eta}\) is a latent random vector of dimension \(n\) which cannot be analytically integrated out of the model.

A common choice of \(\mathbf{R} \left( \boldsymbol{\theta} \right)\) is the Matérn correlation function given by

\[ \begin{align*} R( \mathbf{s}_i, \mathbf{s}_j | \boldsymbol{\theta} ) & = \tau^2 \frac {2^{1-\nu }}{\Gamma (\nu )} \left( {\sqrt {2\nu }} \frac{\|\mathbf{s}_i - \mathbf{s}_j \|}{\phi} \right)^{\nu }K_{\nu } \left( {\sqrt {2\nu }} \frac{\|\mathbf{s}_i - \mathbf{s}_j \|}{\phi} \right). \end{align*} \]

• \(\boldsymbol{\theta} = \left( \nu, \phi \right)'\) where \(\nu\) controls the smoothness of the process and \(\phi\) controls the spatial range of the correlation.

\[ \begin{align*} y(\mathbf{s}_i) | \tilde{\lambda}(\mathbf{s}_i) & \stackrel{iid}{\sim} \operatorname{Poisson}\left( \tilde{\lambda} \left( \mathbf{s}_i \right) \right), \end{align*} \]

where

\[ \begin{align*} \lambda(\mathbf{s}_i) & \sim \operatorname{Gamma} \left( \epsilon_\lambda, \epsilon_\lambda \right). \end{align*} \]

• For the second stage estimates to have finite variance, the first stage posterior distribution \(\left[ \tilde{\lambda} \left( \mathbf{s}_i \right) | \cdot \right]\) should have heavier tails than the target posterior distribution \(\left[ \lambda \left( \mathbf{s}_i \right) | \cdot \right]\).

• \(\tilde{\lambda} \left( \mathbf{s}_i \right) \sim \operatorname{Gamma}(\epsilon_\lambda, \epsilon_\lambda)\) for some small \(\epsilon_\lambda\) will produce a diffuse, analytic first-stage posterior.

  • The analytic solution is beneficial because it allows for fast, efficient simulation of the first stage posterior.

• Metropolis-Hastings update for \(\lambda(\mathbf{s}_i)\) where the proposal distribution \(\lambda^{\star}(\mathbf{s}_i) \sim [\lambda(\mathbf{s}_i) | \mathbf{y}]\) is an independent, discrete sample from the first stage posterior with acceptance probability \(\min(1, a_i^{\star})\) where

\[ \begin{align} \label{eq:MH-poisson} a_i^{\star} & = \frac{[y(\mathbf{s}_i) | \lambda^{\star}(\mathbf{s}_i)] [\lambda^{\star} ( \mathbf{s}_i) | \lambda(\mathbf{s}_{-i}), \boldsymbol{\eta}^{(k_2-1)}, \boldsymbol{\beta}^{(k_2-1)}] [\lambda^{(k_2-1)}(\mathbf{s}_i) | y_i] }{[y(\mathbf{s}_i) | \lambda^{(k_2-1)}(\mathbf{s}_i)] [\lambda^{\star} ( \mathbf{s}_i) | \lambda(\mathbf{s}_{-i}), \boldsymbol{\eta}^{(k_2-1)}, \boldsymbol{\theta}^{(k_2-1)}] [\lambda^{\star}(\mathbf{s}_i) | y(\mathbf{s}_i)] } \\ \nonumber & = \frac{ \left[\lambda^{\star} ( \mathbf{s}_i) | \lambda(\mathbf{s}_{-i}), \boldsymbol{\eta}^{(k_2-1)}, \boldsymbol{\beta}^{(k_2-1)} \right] \left[ \tilde{\lambda}^{(k_2-1)} \left( \mathbf{s}_i \right) \right] }{\left[\lambda^{(k_2-1)} ( \mathbf{s}_i) | \lambda(\mathbf{s}_{-i}), \boldsymbol{\eta}^{(k_2-1)}, \boldsymbol{\beta}^{(k_2-1)} \right] \left[ \tilde{\lambda}^{\star} \left( \mathbf{s}_i \right) \right]}. \end{align} \]

• Requires samples of the latent \(n\)-dimensional vector \(\boldsymbol{\eta}\).

However, the latent parameter \(\boldsymbol{\eta}\) in

\[ \begin{align} \left[\lambda^{\star} ( \mathbf{s}_i) | \lambda(\mathbf{s}_{-i}), \boldsymbol{\eta}^{(k_2-1)}, \boldsymbol{\beta}^{(k_2-1)} \right], \end{align} \]

can be integrated out of the model so that

\[ \begin{align*} \log \left( \lambda(\mathbf{s}_i) \right) | \log \left( \boldsymbol{\lambda}(\mathbf{s}_{-i}) \right) & \sim \operatorname{N} \left( \bar{\mu}_i, \bar{\sigma}^2_{i} \right), \end{align*} \]

where \(\boldsymbol{\lambda}(\mathbf{s}_{-i})\) is the vector of all of the \(\lambda(\mathbf{s}_j)\)s with \(j \neq i\) and

\[ \begin{align*} \bar{\mu}_i & = \mathbf{X}_i \boldsymbol{\beta} + \boldsymbol{\Sigma}_{i, -i} \boldsymbol{\Sigma}^{-1}_{-i, -i} \left( \log \left( \boldsymbol{\lambda} (\mathbf{s}_{-i}) \right) - \mathbf{X}_i \beta \right) \mbox{ and} \\ \bar{\sigma}^2_i & = \boldsymbol{\Sigma}_{i, i} - \boldsymbol{\Sigma}_{i, -i} \boldsymbol{\Sigma}^{-1}_{-i, -i} \boldsymbol{\Sigma}_{-i, i}, \end{align*} \]

with \(\boldsymbol{\Sigma}_{i, i} = \operatorname{Cov} \left( \log \left( \lambda(\mathbf{s}_i) \right), \log \left( \lambda(\mathbf{s}_i) \right) \right)\) and \(\boldsymbol{\Sigma}_{i, -i} = \operatorname{Cov} \left( \log \left( \lambda(\mathbf{s}_i) \right), \log \left( \boldsymbol{\lambda}(\mathbf{s}_{-i}) \right) \right)\).

• Two sample sizes for proof of concept.

• Full rank spatial process of 150 (small) and 300 (medium) samples.

  • Note: typically full-rank spatial models are computationally intractable at about 1000 samples.

  • Computational complexity scales with \(O(n^3)\) and memory scales as \(O(n^2)\).

  • For real-world problems, use reduced rank ideas like those presented earlier to scale to massive datasets.

• Hybrid sampler between MH random sample and IS.
  • Similar to Firefly MCMC.

• Density smooth on first stage samples.
  • Allow for a tuneable random walk rather than independent first stage samples.

• Define transient first stage prior using data.
  • Seems like it would violate rules for specifying the prior.
  • Because transient prior is canceled out in second stage, this issue is irrelevant.

• Coming soon:

  ---

  • BayesMRA – Bayesian spatial and spatio-temporal MRA models (inspired by LatticeKrig).

  ---

  • pgR – Pólya-gamma multinomial regression models.

  ---

  • pollenReconstruction - reconstruction of temperature from pollen.

• Talk available at www.jtipton25.github.io/Los-alamos-2021.

• Website: www.johnrtipton.com.